###################
## JJN Analysis  ##
## 5/29/25       ##
###################
rm(list = ls())

#Loading Packages
P<-c("cjoint", "haven", "ggplot2", "ggtext", "dplyr",
     "ggeffects", "marginaleffects", "AlgDesign", "sandwich", "lmtest", "emmeans", "stargazer")

for (i in 1:length(P)) {
  ifelse(!require(P[i],character.only=TRUE),install.packages(P[i]),
         print(":)"))
  library(P[i],character.only=TRUE)
}

rm(P,i)

#Setting wd
setwd("/Users/morrganherlihy/Library/CloudStorage/Dropbox/Interest_Group_Endorsements/Data_and_Code")
respondent_level_df <- read.csv("respondent_level_vignette_5_29_25.csv")
cj_data_df <- read.csv("respondent_level_conjoint_5_29_25.csv")

#######################
## Study 1: Vignette ##
#######################

#Factoring Vignette Variables for Analysis
respondent_level_df$Treatment <- factor(respondent_level_df$Treatment, 
                                        levels = c("Forthcoming", "Avoid", "Client"))
table(respondent_level_df$Treatment)

respondent_level_df$PID <- factor(respondent_level_df$PID, 
                                  levels = c("Democrat", "Independent", "Republican"))
table(respondent_level_df$PID)

respondent_level_df$Hispanic <- factor(respondent_level_df$Hispanic, 
                                       levels = c("non-Hispanic", "Hispanic"))
table(respondent_level_df$Hispanic)

respondent_level_df$Female <- factor(respondent_level_df$Female, 
                                     levels = c("Male", "Female"))
table(respondent_level_df$Female)

respondent_level_df$PolInt <- factor(respondent_level_df$PolInt, 
                                     levels = c("Not Interested", "Interested"))
table(respondent_level_df$PolInt)

respondent_level_df$Religious <- factor(respondent_level_df$Religious, 
                                        levels = c("Not Religious", "Religious"))
table(respondent_level_df$Religious)

#Model Results, Base Model
m1_vignette <- lm(Nom_SupportDV ~ Treatment + Prochoice_Op, respondent_level_df)
summary(m1_vignette)

##Model Results, Interaction Model
m2_vignette <- lm(Nom_SupportDV ~ Treatment * Prochoice_Op, respondent_level_df)
summary(m2_vignette)

##Model Results, Full Model
m3_vignette <- lm(Nom_SupportDV ~ Treatment * Prochoice_Op + PID + Hispanic + Female + Age + PolInt + Ideology + Religious, respondent_level_df)
summary(m3_vignette)

##Table 1, Regression Summary from Vignette Models
models <- list(m1_vignette, m2_vignette, m3_vignette)

stargazer(models,
          header = FALSE,
          type = "latex",
          title = "Survey Experiment Conditions",
          dep.var.labels.include = FALSE,
          dep.var.caption = "",
          model.names = FALSE,
          model.numbers = FALSE,
          column.labels = c("Base", "Interaction", "Full"),
          omit.stat = c("LL","ser","f"),
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes.append = FALSE,
          out = "/Users/morrganherlihy/Desktop/Judging_Judicial_Nominees/figures/PRQ/table1_main.tex")

##Figure 1, Predicted Support across Treatments, Vignette
#Predict
predictions_interaction <- predictions(m2_vignette, newdata = datagrid(
  Treatment = unique(respondent_level_df$Treatment),
  Prochoice_Op = seq(min(respondent_level_df$Prochoice_Op, na.rm = TRUE),
                     max(respondent_level_df$Prochoice_Op, na.rm = TRUE),
                     length.out = 100)))

#Set Plot Theme, Vignette
jjn_theme <- theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
        plot.title = element_text(hjust = 0, size = 8, face = "bold"), 
        axis.title = element_text(hjust = 0.5, size = 8, face = "bold", margin = margin(r = 10)), 
        axis.text = element_text(size = 7), 
        axis.ticks = element_line(size = 0.5, color = "black"), 
        plot.caption = element_markdown(hjust = 0, size = 8),
        plot.margin = margin(t = 0, b = 0, l = 10, r = 10, unit = "pt"), 
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7),
        legend.position = "bottom",               
        legend.direction = "horizontal",       
        legend.box.just = "left", 
        legend.key.width = unit(1.75, 'cm'),
        legend.key.height = unit(.5, "cm"),
        legend.margin = margin(t = 5, r = 0, b = 5, l = 0), 
        legend.justification = c(0, 1))

setwd("/Users/morrganherlihy/Desktop/Judging_Judicial_Nominees/figures/PRQ")
pdf("vignette_interaction_5.25.25.pdf",width=5,height=3.5,paper='special')
ggplot(predictions_interaction, aes(x = Prochoice_Op, y = estimate, linetype = factor(Treatment))) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(Treatment)), alpha = 0.3) +  
  labs(x = "",
       y = "Predicted Support for Confirmation",
       linetype = "Treatment",
       fill = "Treatment") + 
  geom_line(size = 1) +  
  scale_linetype_manual(values = c("dotted", "solid", "dashed")) +  
  scale_fill_manual(values = c("gray30", "gray50", "gray70")) + 
  scale_x_continuous(
    limits = range(predictions_interaction$Prochoice_Op),
    labels = c("Strongly\nPro-Life", "Somewhat\nPro-Life", 
               "Somewhat\nPro-Choice", "Strongly\nPro-Choice")) +  
  scale_y_continuous(limits = c(1.75, 3.75), breaks = c(2, 2.5, 3, 3.5)) +  
  jjn_theme
dev.off()

#Differences between Treatments for In-Text Comparisons
compare_estimates_across_prochoice <- function(data, prochoice_values, ref_group, compare_groups) {
  results <- data.frame(ProChoice_Op = numeric(), Comparison = character(), 
                        Difference = numeric(), SE_Diff = numeric(), 
                        t_value = numeric(), p_value = numeric(), 
                        stringsAsFactors = FALSE)
  
  for (prochoice_value in prochoice_values) {
    df_filtered <- data %>% filter(Prochoice_Op == prochoice_value)
    
    estimates <- df_filtered %>% 
      dplyr::select(Treatment, estimate, conf.low, conf.high)
    
    estimates <- estimates %>%
      mutate(SE = (conf.high - conf.low) / (2 * 1.96))
    
    ref_est <- estimates %>% filter(Treatment == ref_group)
    ref_estimate <- ref_est$estimate
    ref_SE <- ref_est$SE
    
    for (group in compare_groups) {
      comp_est <- estimates %>% 
        filter(Treatment == group)
      
      diff_estimate <- ref_estimate - comp_est$estimate
      SE_diff <- sqrt(ref_SE^2 + comp_est$SE^2)
      
      t_value <- diff_estimate / SE_diff
      p_value <- 2 * pt(-abs(t_value), df = Inf)  
      
      results <- rbind(results, data.frame(ProChoice_Op = prochoice_value,
                                           Comparison = paste(ref_group, "vs", group),
                                           Difference = diff_estimate,
                                           SE_Diff = SE_diff,
                                           t_value = t_value,
                                           p_value = p_value))
    }}
  return(results)
}

compare_results_forthcoming <- compare_estimates_across_prochoice(
  data = predictions_interaction, 
  prochoice_values = c(1, 2, 3, 4),  
  ref_group = "Forthcoming", 
  compare_groups = c("Client", "Avoid"))

compare_results_avoid <- compare_estimates_across_prochoice(
  data = predictions_interaction, 
  prochoice_values = c(1, 2, 3, 4),  
  ref_group = "Avoid", 
  compare_groups = c("Client", "Forthcoming"))

compare_results_client <- compare_estimates_across_prochoice(
  data = predictions_interaction, 
  prochoice_values = c(1, 2, 3, 4),  
  ref_group = "Client", 
  compare_groups = c("Forthcoming", "Avoid"))


#######################
## Study 2: Conjoint ##
#######################

#Factoring Conjoint Variables for Analysis
cj_data_df$Abortion <- factor(cj_data_df$Abortion, levels = c("Pro-Life", "Pro-Choice"))
table(cj_data_df$Abortion)

cj_data_df$A_Gender <- factor(cj_data_df$A_Gender, levels = c("Male", "Female"))
table(cj_data_df$A_Gender)

cj_data_df$B_Law.School <- factor(cj_data_df$B_Law.School, levels = c("Top-Tier", "Mid-Tier", "Lower-Tier"))
table(cj_data_df$B_Law.School)

cj_data_df$C_Occupation <- factor(cj_data_df$C_Occupation, levels = c("Federal Judge", "State Supreme Court Justice", "Assistant U.S. Attorney", "Public Defender", "Private Attorney", "ACLU Attorney", "Law Professor"))
table(cj_data_df$C_Occupation)

cj_data_df$D_IG.Statement1 <- factor(cj_data_df$D_IG.Statement1, levels = c("Neutral","Support","Oppose","Qualified"))
table(cj_data_df$D_IG.Statement1)

cj_data_df$E_Nominee.Quote1 <- factor(cj_data_df$E_Nominee.Quote1, levels = c("Evasive", "Forthcoming"))
table(cj_data_df$E_Nominee.Quote1)

#Theme for Conjoint Plots: 
conjoint_theme <- theme(
  panel.border = element_rect(colour = "black", fill = NA, size = 1.25),
  panel.background = element_rect(fill = "white", colour = "white", size = 0.5),
  panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid', colour = "grey80"),
  legend.position = "none",
  strip.background = element_rect(colour = "black", fill = NA, size = 1.25),
  strip.text.x = element_text(size = 12, color = "black", face = "bold", margin = margin(t = 4, b = 4)),
  axis.text.y = element_text(hjust = 0, size = 11),
  axis.text.x = element_text(colour = "black", size = 11),
  axis.title = element_text(size = 14),
  plot.margin = margin(t = 15, r = 15, b = 15, l = 15))

#Conjoint, Base Results 
m1_cj <- amce(selected ~ A_Gender + B_Law.School + C_Occupation + D_IG.Statement1 + E_Nominee.Quote1, 
              data = cj_data_df, respondent.id = "Case_ID")
summary(m1_cj)

#Figure 2, Conjoint Base Results and Figure 2
setwd("/Users/morrganherlihy/Desktop/Judging_Judicial_Nominees/figures/PRQ")
pdf("cj_base_5.25.25.pdf", width = 7, height = 7,paper='special')
plot(m1_cj, plot.theme=conjoint_theme, attribute.names = c("Gender", "Law School Ranking", "Occuption", "Interest Group Statement", "Nominee Quote"), xlab = "Change in Likelihood of Selection")
dev.off()

#Conjoint, Conditional Results and Figure 3
m2_cj <- amce(selected ~ A_Gender * Abortion + B_Law.School * Abortion + C_Occupation * Abortion + D_IG.Statement1 * Abortion + E_Nominee.Quote1 * Abortion, data = cj_data_df, respondent.id = "Case_ID", respondent.varying = "Abortion")
summary(m2_cj)

pdf("cj_cond_5.25.25.pdf",width=7,height=7,paper='special')
plot(m2_cj, plot.theme = conjoint_theme, plot.display = "interaction", attribute.names = c("Gender", "Law School Ranking", "Occuption", "Interest Group Statement", "Nominee Quote"), xlab = "Change in Likelihood of Selection")
dev.off()

#Conjoint, Conditional and Interaction Results and Figure 4
m3_cj <- lm(selected ~ A_Gender * Abortion + B_Law.School * Abortion + C_Occupation * Abortion + E_Nominee.Quote1 * Abortion * D_IG.Statement1, data = cj_data_df)
coeftest(m3_cj, vcov = vcovCL(m3_cj, type = "HC1",cluster = cj_data_df$Case_ID))

pred_m3_cj <- ggpredict(model=m3_cj,
                        terms=c("D_IG.Statement1 [all]", "E_Nominee.Quote1 [all]", "Abortion [all]"), 
                        vcov.fun="vcovCL",
                        vcov.type="HC1",
                        vcov.args=list(cluster=cj_data_df$Case_ID))    

pred_m3_cj <- pred_m3_cj %>% 
  rename(D_IG.Statement1 = x,
         Predicted = predicted,
         LowerCI = conf.low,
         UpperCI = conf.high,
         E_Nominee.Quote1 = group,
         Abortion = facet) 

#Plot new predictions with clustered standard errors
jjn_theme <- theme_minimal() +
  theme(panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
        plot.title = element_text(hjust = 0, size = 12, face = "bold"), 
        axis.title = element_text(hjust = 0.5, size = 12, face = "bold", margin = margin(t = 10, r = 10, b = 10, l = 10)), 
        axis.title.x = element_text(hjust = 0.5, size = 12, face = "bold", margin = margin(t = 15)),
        axis.title.y = element_text(hjust = 0.5, size = 12, face = "bold", margin = margin(r = 15)),
        axis.text = element_text(size = 8), 
        axis.ticks = element_line(linewidth = 0.5, color = "black"), 
        plot.caption = element_markdown(hjust = 0, size = 12),
        plot.caption.position = "plot", 
        plot.margin = margin(t = 0, b = 0, l = 10, r = 10, unit = "pt"), 
        strip.text = element_text(size = 14, face = "bold"),
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.position = "bottom",               
        legend.direction = "horizontal",       
        legend.box.just = "left", 
        legend.margin = margin(t = 5, r = 0, b = 5, l = 0), 
        legend.justification = c(0, 1))

pdf("cj_cond_interaction_5.25.25.pdf",width=5,height=3.5,paper='special')
ggplot(pred_m3_cj, aes(x = D_IG.Statement1, y = Predicted, group = E_Nominee.Quote1, shape = E_Nominee.Quote1, 
                       color = E_Nominee.Quote1)) +
  geom_point(position = position_dodge(width = 0.3), size = 3) +  
  geom_errorbar(aes(ymin = LowerCI, ymax = UpperCI), 
                position = position_dodge(width = 0.3), width = 0) + 
  facet_wrap(~Abortion) + 
  scale_shape_manual(values = c("Forthcoming" = 19, "Evasive" = 17)) +  
  scale_color_manual(values = c("Forthcoming" = "gray40", "Evasive" = "gray60")) +
  labs(title = "",
       x = "Interest Group Statement",
       y = "Predicted Probability",
       shape = "Nominee Quote", 
       color = "Nominee Quote") +
  jjn_theme
dev.off()

################
## Appendices ##
################

#Appendix B, Direct Effects Model
##Model, Base
m1_vignette_appendix <- lm(Nom_SupportDV ~ factor(Treatment), respondent_level_df)
summary(m1_vignette_appendix)

##Model, Full Model
m2_vignette_appendix <- lm(Nom_SupportDV ~ factor(Treatment) + PID + Hispanic + Female + Age + PolInt + Ideology + Religious, respondent_level_df)
summary(m2_vignette_appendix)

##Appendix B Table, Regression Summary from Vignette Models
models_appendix <- list(m1_vignette_appendix, m2_vignette_appendix)

stargazer(models_appendix,
          header = FALSE,
          type = "latex",
          title = "Survey Experiment Conditions",
          dep.var.labels.include = FALSE,
          dep.var.caption = "",
          model.names = FALSE,
          model.numbers = FALSE,
          column.labels = c("Base", "Interaction", "Full"),
          omit.stat = c("LL","ser","f"),
          digits = 2,
          digits.extra = 2,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes.append = FALSE, 
          out = "/Users/morrganherlihy/Desktop/Judging_Judicial_Nominees/figures/PRQ/tableB1_appendix.tex")

#Appendix C, Conjoint Comparing Marginal Means Tables C1 and C2
clustered_vcov <- vcovCL(m3_cj, 
                         cluster = cj_data_df$respondent, 
                         type = "HC1") 
marg_means <- emmeans(m3_cj, ~ E_Nominee.Quote1 * Abortion * D_IG.Statement1, 
                      vcov. = clustered_vcov)
summary(marg_means)

#Appendix C, Conjoint Comparing Marginal Means p-values, Tables C3 - C6
contrast(marg_means, method = "pairwise", adjust = "none")
cont_mm <- summary(contrast(marg_means, method = "pairwise", adjust = "none"))

